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NUMERICAL PREDICTION OF TURBULENT THREE-DIMENSIONAL 


JUNCTURE REGION FLOW USING THE COMOC III 
3DPNS COMPUTER PROGRAM 

By 

A. J, Baker, P, D. Manhardt, 

And J, A. Orzechowski 
Computational Mechanics Consultants, Inc. 
Knoxville, TN 37920 


SUMMARY 


A numerical solution algorithm is established for prediction of subsonic 
turbulent three-dimensional flows in aerodynamic configuration juncture 
regions. In concert with a full three-dimensional exterior potential flow 
solution, the developed parabolic algorithm yields prediction of the details 
of the corner region flows. A turbulence closure model is established using 
the complete Reynolds stress. Pressure coupling is accomplished using the 
concepts of complementary and particular solutions to a Poisson equation. 
Numerical results are presented for prediction of the three-dimensional 
turbulent flow in the juncture of two intersecting parabolic arc airfoils. 
Specifications for data input juncture geometry modifications are presented. 
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INTRODUCTION 


This effort is in support of the advanced aerodynamics area of the Air- 
craft Energy Efficiency (ACEE) program. Specifically, a viscous-inviscid flow 
analysis is required developed and evaluated for prediction of the viscous and 
turbulence effects in a geometry corresponding to juncture regions of wing- 
body, pylon-wing, and/or winglet-wing intersections. The effort reported 
herein focuses on the development and formulation of a method and an associated 
computer program for calculation of the viscous and turbulent flow in the 
region of wing-body, pylon-wing, and winglet-wing intersections considering 
the geometries appropriate to supercritical wing design (including winglets), 
and the interference effects due to propulsion installation. 

The state-of-the-art of three-dimensional boundary-layer calculations has 
matured to where finite wings as well as bodies can be treated. Prior to the 
current analysis, no valid technique has been developed which can properly 
analyze the juncture region of a wing-body combination. This void in analysis 
capabilities handicaps the designing of aircraft free from flow separation and 
the evaluating and minimizing of aircraft drag. As a consequence, current 
three-dimensional boundary- layer programs require inboard starting profiles 
that must either be obtained from an infinite swept-wing solution or an 
infinite tapered wing solution, otherwise an assumption of symmetry at the 
starting line is required. Similar assumptions are required at the tip. 

Three-dimensional flows over finite intersecting surfaces fall into two 
general categories as a function of the flow structure in the plane transverse 
to its predominant direction. This in turn is strongly dependent on the flow 
Reynolds number. The essential geometry of the juncture problem corresponds 
to flow in the immediate vicinity of one corner of a duct, see Figure 1. 

The solution domain is assumed unbounded in the first quadrant of the transverse 
plane. Carrier (ref. 1) formulated the incompressible laminar flow 
problem in boundary layer similarity form, and expressed the corner effects 
as an alteration to the Blasius solution. This analysis was incorrect however, 
since it failed to account for transverse plane vorticity. Rubin (ref. 2) 
formulated the low-speed laminar flow corner problem in completeness, using 
resolution of the flow domain into potential, boundary layer and corner layer 
regions. Subsequently, Pal and Rubin (ref. 3) and Rubin and Grossman (ref. 4) 
evaluated asymptotic characteristics and numerical solutions for incompressible 
laminar flow. Extension to laminar compressible flow is reported by Weinberg 
and Rubin (ref. 5). Ghia and Davis (ref. 6, 7) and Ghia (ref. 8) evaluated 
use of optimal coordinates for the compressible laminar flow corner. In all 
cited instances, numerical solutions of the governing equations were 
established for a symmetric half domain formed by the bisector of the right 
angle corner, one wall and asymptotic approach to two-dimensionality. 
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Tokuda (ref. 9) observed that the isovels computed by Rubin and Grossman 
are incorrect in comparison to the experimental data of Zamir and Young 
(ref. 10). This was attributed to lack of account of a fourth region tucked 
inside the corner layer, and he performed a Stoke 's flow expansion therein 
and theoretically combined solutions by overlapping various regions. However, 
the experimental data is for turbulent flow, and Zamir and Young note a funda- 
mental difference in the velocity contour bulges nearby the corner. As shown 
in Figure 2, laminar corner flow apparently induces mass entrainment directed 
along the walls into the corner with efflux along the bisector. For turbulent 
flow, the data indicates influx into the corner along the bisector with out- 
ward transport along the walls. Hence, the fundamental mechanisms for laminar 
and turbulent corner interaction must be distinct, and analysis of the latter 
is the present primary requirement. 

Bragg (ref 11) analyzed turbulent incompressible corner flow and deter- 
mined the distribution of Turbulence Kinetic Energy, TKE, in the corner. He 
observed that symmetry in a turbulent flow can occur at best only under the 
most controlled laboratory conditions, indicating that a full domain numerical 
solution is required. Eichel brenner and Preston (ref. 12) present a compre- 
hensive theoretical analysis of the role of secondary flow (in the transverse 
plane) on turbulent corner flow. They propose existence of systems of secon- 
dary vortices on the plate, see Figure 3, with the resultant corresponding 
pressure distribution linked to deformation in the normal components of the 
Reynold's stress tensor. They conclude that the resultant anisotropy of the 
stress is the source of the observed secondary flow behavior that induces mass 
flux into the corner along the bisector, see Figure 4. They further note the 
resultant vorticity intensity is highest if allowed to develop freely as in 
the unbounded juncture region. Gessner (ref. 13) refutes this analysis, pre- 
senting data which fails to indicate the required transverse pressure undula- 
tion. He presents a detailed theoretical analysis of mechanisms inducing 
secondary flow; an energy balance of the mean flow isolates two dominant terms 
responsible for observed corner influx along the bisector. Experimental data 
indicate that the gradient, normal to the corner bisector, of the corresponding 
Reynolds shear stress component, ufun , is the dominant mechanism inducing and 
maintaining the corner flow. The induced secondary mean flow convects axial 
mean velocity, vorticity and mean flow energy into the corner as a first-order 
effect. The axial vorticity component is a second-order effect as is the 
energy of the turbulence velocity field and its convection by the secondary 
mean flow. 

The essential flow character in the juncture region thus appears the 
result of a delicate balance between turbulence phenomena and the induced 
secondary mean flow velocity field. Anisotropy of the Reynold's stress tensor 
is important, but apparently not a dominant influence, but the distribution 
of the shear components in the immediate corner region appears fundamentally 
important. With this insight, a parabolized form of the three-dimensional, 
time-averaged Navier-Stokes equations for prediction of steady, turbulent, 
compressible shock-free flow in a juncture region is derived. Determination 
of the six components of the Reynolds stress field is accomplished using a 
parabolic simplification to a constitutive equation. The resultant tensor 
field is anisotropic and requires solution of parabolized three-dimensional 
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forms of the transport equation for turbulence kinetic energy and isotropic 
dissipation function. 

The developed viscous and turbulent flow description forms one part of a 
viscous-inviscid interaction algorithm for determination of the complete 
juncture region flowfield. An assumed "viscous displacement distribution" is 
used to generate an initial surface estimate for a three-dimensional potential 
flow solution. This yields the corresponding surface pressure distribution on 
the freestream closure segment of the viscous flow solution domain. Pressure 
coupling between the two solutions is accomplished through the complementary 
solution to a Poisson equation. The pressure particular solutions yields the 
detailed distribution near the juncture waVTs in concert with computational 
enforcement of conservation of mass. Following the initial solution, the 
computed viscous flow freestream velocities can be employed as onset velocities 
for a subsequent three-dimensional potential flow refinement, and an iterative 
interaction algorithm thereby developed. 

A computer program has been written for solution of the developed three- 
dimensional Navier-Stokes (3DPNS) equation system. The Hess (ref. 14) computer 
program was employed to generate the inviscid flow pressure distribution on the 
juncture region surface of two intersecting non-lifting 10% thick parabolic 
arc airfoils with coincident leading edges and zero sweep. The 3DPNS system 
was solved and the three-dimensional distributions of mean velocities, Reynolds 
stress components, pressure distributions and onset flow refinements for the 
potential flow refinement determined. The results of this study are reported 
herein, as well as input instructions for use of the developed 3DPNS option of 
the COMOC III compute program for juncture region flow prediction. 


SYMBOLS 

a boundary condition coefficient; direction cosines 

A Van Driest damping function; constant 

b body force 

c isentropic sound speed; coefficient 

C coefficient 

Cp pressure coefficient 

C skin friction coefficient 

P 

d differential 

e fintie element index 


5 


f 


f 

Fr 

g 

h 
H 
k 
k 


A- 


function of known argument; coordinate c'H'jVe 
Froude Number , 

gravity acceleratioh; function; source term 
metric coefficient; mesh parameter; integration step size 
boundary layer shape factor 
thermal conductivity; turbulence Up^tic energy 
generalized diffusion cpefficf"^ 
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. V . . Re 
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Xi 


differential operator; lengt^^|®Te'''^‘ '\- 
differential operator; length 
Mach number; number of finite elements spanning 
unit normal vector; nodes per element; dimensionality 
finite element interpolation function 
pressure; generalized parameter; iteration index 
Prandtl number 

generalized dependent variable 
generalized discretized dependent variable 
Ifbmain of elliptic operator 
Reynolds number 

finite element assembly operator; surface 
time 




■■ 






V*v* 


velocity vector, 
reference velocity 
friction velocity 
computational velocity vector 
Cartesian coordinate system 
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Y ratio of specific heats 

3 R closure of solution domain R 

6 Kronecker delta, boundary layer thickness; increment 

6 * boundary layer displacement thickness 

A increment; element measure 

e turbulence dissipation function 

0 boundary layer momentum thickness 

transformed coordinate system 
transformed coordinate system 
K Karman coefficient (MLT) 

A multiplier; turbulence sublayer constant 

y dynamic viscosity 

V kinematic viscosity; general diffusion coefficient 

p density 

a- ■ mean flow Stokes stress tensor 

* %J 

T Reynolds stress tensor; wall shear; integration kernel 

<fi Velocity potential function 

X generalized initial -val ue operator 

finite element natural coordiante system 
0) turbulence damping factor 

global solution domain 

Superscripts: 

e effective value 

t dimension of R 


t 


turbulent 



T matrix transpose 

+ turbulence correlation function 

mass-weighted time-average 

time-averaged 

unit vector 

^ mass-weighted fluctuating component; ordinary derivative 

Subscripts: 

“ global reference condition 

e finite element domain 

i,j»k,£ tensor indices 

non-tensor index 

I freestream reference condition 

n normal 

0 initial condition; stagnation reference 

t time derivative 

w wall reference condition 

Notation: 

{ } column matrix 

[ ] square matrix 

(J union 

A intersection 

E belongs to 


8 


PROBLEM DESCRIPTION 


Governing Differential Equation System 

The description of a state point in multi -dimensional fluid mechanics is 
contained within solution of a coupled nonlinear partial differential equation 
system describing conservation of mass, momentum and energy. Unique solutions 
are obtained upon closure by specification of constitutive behavior and 
boundary conditions. In Cartesian tensor notation, the non-dimensional 
conservation form of the Navier-Stokes equation system is 



= If 


9 

dx 


— [pu-u. + p6.. - Re’^a. .1 + Fr"^pb. 


’’‘JL 


pUjH - (y - - Re“^Pr'^ 


9H 

dx 




( 1 ) 

( 2 ) 

(3) 


The dependent variable in equations (l)-(3) have their usual interpreta- 
tion in fluid mechanics, i.e., p is mass density, U4 is the velocity vector, 
p is the static pressure, b is a body force, and H ^ is the stagnation 
enthalpy. Furthermore, ajj is the Stokes stress tensor, defined as 








3 3X|j ”ij 


( 4 ) 


where y is the dynamic viscosity. The equation of state for a perfect gas 
closes the definition, and the non-dimensional parameters of impact are: 


Pc^U^£ 

Reynolds Number: Re h 


(5) 
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Prandtl Number; 


Pr i 


( 6 ) 



Mach Number: M e — 

oo c 


(7) 


Herein, £ is a characteristic scale length and c is the isentropic sound 
speed. 

Equations (l)-(7) are valid for both laminar and turbulent flows. For 
the latter, however, their solution becomes tractible in a practical sense 
only after time-averaging, and mass-weighted time-averaging is assumed to 
serve present requirements. Therefore, the Reynolds decomposition is defined 
as (ref. 15) 




Uj(x^,t) = + ur(x.,t) 


(8) 


The mass-weighted time-average velocity is defined as 



This operation yields the important relation 



Equations (8)-(9) are also employed to define the time-averaged and fluctuating 
components of enthalpy as 


H = ^ 2 V 1 




(12) 


Substitution of equations (8)-(9) into (l)-(4), time-averaging and 
collecting terms yields the time-averaged Navier-Stokes equation system 








8(PU.) 

— 5T- * sf: 


+ P6^J + pu,u'J = 0 




im - ^ ^ ^ [=j(3H) * ?iruj.-(Y - 1)«^> (s,5,j . 

, ^ 


(13) 

(14) 


- Re“iPr~^ij 

A. 



(y - 1 )M^ 1^ * 0 


(15) 


The three-dimensional parabolic Navier-Stokes (3DPNS) equations are 
required established to describe the steady time-averaged viscous and turbulent 
flow of a compressible fluid in the juncture region. The assumptions for 3DPNS 
are: 


(1) A predominant mean flow direction is uniformly discernible. 

(2) In this direction (only) diffusion processes are negligible 
compared to convection, and, 

(3) Overall three-dimensional elliptic character is provided by 
interaction with the potential freestream flow. 


The parabolic approximation to equations (13)-(15) basically constrains 
summation limits. Assuming the Xj coordinate aligned with the direction of 
predominant flow, the parabolic approximation is concisely expressed as 




3u^ 

J 


3u, 


3x. 


2.!^ 

3 3x,. 




( 16 ) 
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The subscript bar notation denotes the index not eligible for the 
summation convention, but merely takes on the value of the synonymous tensor 
index. The parabolic approxima^’i'''' equation (15) also requires eliminating 
the Xj comoonon^- ■■ i involving the time-averaged viscosity 

IT . The I ^ 1 ~ u + ^ value character on the Xj 

coordinan ^ s h + 2 V ^ ields the desired form upon ide ntifi - 

cation of| ..tc moaei for the Reynolds stress tensor - pu:[Uj . 


Constitutive Closure For Turbulence 


The 3DPNS equation system becomes closed upon identification of the 
components of the Reynolds stress tensor. Based upon concepts in continuum 
mechanics, c.f. Lumley (ref, 16), a constitutive equation for the kinematic 

Reynolds stress tensor - utu': is in the form 

* J 


- u.u: 
1 0 





C ^ 
e 




^ + 







(17) 


The coefficients result from if'eexpression of triple correlations within the 
Reynold's stress transport equation using the model of Launder, Reece and 
Rodi (ref. 17). It is a generalization of the original analysis by Gessner 
and Emery (ref. 18) for an equilibrium steady three-dimensional parabolic 
flow. In equation (17) a-jj is a diagonal tensor defined as 




(18) 


Hence, equation (17) represents an anisotropic tensor as indicated required 
in the analysis of reference 12. In equation (18) the a-,- are the coefficients 
that admit anisotropy and ax = Cx>a 2 = a 3 = C 3 is suggested, where the C-j are 
defined by Launder, Reese and Rodi (ref. 17)as 


i 

I 

I 

! 


krl- 22(C,1 - 1) - 6(4C,; - 5) 
33(C„ - 2C,J 1 


- 2C,2) 




C 3 ^ 


22(Cji - 1) - 12(3C|2 - 1-) 
33(C^i - 2C,2} 


' * 


^ _ 44Cei - 22C<iiCi2 - 128C02 - 36Cj2 + 10 . 
^ " 165(Cpi - 26^2)^ S 


(19) 
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In equations (19) C/Si and Zf,x are "universal" empirical- constants derived 
by Hanjalic and Launder (ref. 19); suggested values are C)ji = 2.8 and 
Cjii 2 « 0.45. Additionally, in equation (17). k is the kinetic energy of the 
turbulence velocity fieiH 


H = h + ? 2 ’ 

' 1 . 


(20) 


and e is tbe isotropic dissipation function defined by the contraction 


9ur9u: 

3 


k"''k 

The governing differential equations for k and e are (ref. 15) 


(21) 




le 


e e "i jm 


'2 C^ . 
'e k 


The various C|^ and C“ are empirical model constants to be specified. 


( 22 ) 


(23) 


Three-Dimensional Parabolic Navier-Stokes Equations 

The 3DPNS equation system for solution of isoenergetic juncture region 
flow is obtained by deletion of the time derivatives in equations' (13)-(14) 
and (22) -(23) and eliminating the Xj stress components in equation (17) in 
the manner illustrated for the Stoke 's stress tensor. The resultant 3DPNS 
equation system requiring numerical solution is 
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(24) 


/ T \ 
1 


= h + ^ UiU. + 2 





„ -^'{pu,u ) + — -B- _ 9 


JL 5-7 

Re ^- P^^il 




= 0 




(25) 


^^s#.;-3-!' ' , , , _ ,, 

L(pU 2) = ^ (pU^Uo) 



1 »? 
I C^\/ 





(26) 




.a *f v li . ^ 

13' 3Xg 




i. fliil 


■R6 3 j 






(27) 



„3x^ 3x^ 


R«v' 


t ^ 


9k 


9x 


3u, 


i.fi'i'Hj;^,.* = i (28) 




'-(e) - 3^ 


.-i'v 

.1 ^ 
-A 





0^'" 


(29) 


Equations (24)-(29) introduce the 3DPNS limited summation index convention 
1 £(i, j)£3, 2 £(k, Ji)£ 3. The dilitation term in the Stoke's stress tensor 
has been deleted from equations (25)-(27) as negligibly small. 

The parabolic approximation to equation (17) for use in completing terms 
in equations (25)-(29) is (ref. 16). 
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For equation (25) the divergence term becomes 




^ '"'1^ 


3x 


^ "3^ 


P^^S, 


. -,9 

»%x , 


(30) 


Re 


-■21 1 




Defining an "effective" kinematic viscosity 




(32) 


equation (31) becomes 


3X( 


- 3u, 

_y_ 1 


Re 3Xji " 


* 

3 

[- eff 5“il 

- 3Xj 

3X 


which corresponds exactly to the familiar boundary layer torm. 


(33) 


In the transverse flow momentum equations (26)-(27) the divergence terms 
become, using equations (30) 
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The first right-hand term is basically identical to equation (33) for Uj , 
with the addition of terms involving shear of the alternative transverse 
plane mean velocity component u . The second term is a non-homogeneity 
involving only shear components^of the predominant velocity ui which acts as 
a source term within the transverse momentum equations. 

The divergence term common to both the k and e equations becomes 


1 9 

ic,p Jiucu: M 

9 

[2 c 5 J^ ial 

■ 

4 . 1- 

3 q e 1 > 

« _ 

X 

3 


(35) 


where q represents either k or e. Equation (35) is identical to the Reynold's 
stress contribution in equation (32) with C 4 replaced by Cq. The production 
terms for both k and e are basically identical, with the latter multiplied by 
e/k. The Reynold's stress-velocity shear contraction under the parabolic 

approximation becomes 






The second term is assumed negligible due to the appearance of 

8Ui/8Xi << 9Ui/9X^. 


Pressure Resolution and Mass Conservation 


A key element for the viscous-inviscid interaction algorithm for juncture 
region flow is the pressure coupling between the'two fidwfield descriptions. 
v^golution of the full three-dimensidhal linear ISplaci an equation governing 
^e inviscid, assumed-irrotational.^ul)Sonic exterior fjto^dan be readily 
accomplished with current computer pr^^Tg^wis , e . Thij;;is possibte^i^!^^^ 
using the computational equivf^dnt of eilher flow tSfT^ehcy boundary conditions 
(on a stream surface) or specification of the "onset" velocity on any surface^.^'^.^^'Y 
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This is defined as the difference between the local velocity vector and 
reference freestream vector. Many two-dimensional interaction algorithms 
assume that the (aerodynamic) surface augmented by the boundary layer displace- 
ment thickness 6* corresponds to the inviscid streamline. However, since 6* 
could be difficult to define for other than an elementary juncture geometry, 
the specifications and iterative refinement of the inviscid onset velocity 
appears more attractive. 

This is the sole boundary condition requirement for determination of the 
three-dimensional inviscid flow pressure distribution (Cp) on the viscous flow 
domain intersection with the freestream. The point of departure for derivation 
of the pressure coupling algorithm is the steady parabolized form mean flow 
momentum equations, (25)-(27). It is based on the two-dimensional transverse 
plane Poisson equation obtained by application of the divergence operator. 
Recalling the limited summation convention, 1 j ji 3, and 2 £ k, £ £ 3, this 
equation takes the form 



The algorithm is based upon the observation that the solution of a linear, 
elliptic Poisson equation consists of a complementary and a particular con- 
tribution as 



(38) 


By definition, the compl ementary solution satisfies the homogeneous part of 
equation (37), i.e., 

^ o • ' > 






(39) 


subject to the known boundary conditions for p(x^). Since the bounding 
inviscid flow pressure distribution is everywhere known. 



(40) 


In equation (40), Xj^ indicates Xj^ constrained to the boundary of the viscous 
solution domain, = 3 R 2 x x^, and pi is the (inviscid flow) pressure level 
distribution. Elsewhere on the appropriate boundary condition is vanish- 
ing normal gradient, i.e.. 


*(Pc> = 






(41) 
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i%f. 


The particular solution Pp is any function satisfying equation (37);'^ith'''' 
homogeneous boundary conditions on closure segments of fi coincident wibh-ijl 
known. Elsewhere, the boundary condition for pp is established from the inner 
product of the parent steady momentum equations (26)-(27) with the closure i 
surface normal. Since the convection term vanishes identically at a waTl 'or 
symmetry plane, the appropriate form is 




.*.5^ 

«(Pp) 


Sfe 




^.n 
9Xj; ^ 


8Xk 






Furthermore, since retaining the second term in equation (42) obliterates th^ 
appearance of the boundary condition stress term in equation (37), the appfo-^’ 
priate Poisson equation for Pp must be the inviscid form, ke., • -i- 




L(Pi 


ax. 


3x, 




= 0 


(43) 




Equations (37)-(43) define the algorithm for determination of the pressure 
distributions within the three-dimensional parabolic viscous flow solution 
domain. The overall elliptic character of the subsonic flowfield is enforced 
by the fully three-dimensional exterior potential flow pressure distribution 
boundary conditions. Hence, a sequential iteration between the two problem 
specifications is required. To avoid a possible destabilizing of the parabolic 
equation solution, the axial pressure gradient distributions driving u^ are 
assumed constituted solely of the complementary pressure field. The current 
computed particular and complementary pressures are summed to form the gradients 
in the plane transverse to Ui. Following the first sweep of the parabolic 
solution, a refined complementary field is constructed from the three-dimen- 
sional potential flow solution on the modified viscous surface. For the second 
parabolic solution, the current complementary and previous-iterate particular 
pressure distribution are added to form the axial pressure gradient for Ui. 
Hence, in convergence of the interaction algorithm, the current-iterate par- 
• ticular pressure distribution should become negligibly small. 

• % 

Equation (43) can be simplified by insertion of the continuity equation 
(24). However, a direct enforcement of conservation of mass is also required 
established. Hence, assume the solution mean flow velocity fieTd u. comprised 
of the sum of the solution v-, of the steady flow momentum equation (25)-(27), 
and an irrotational correction velocity field determined frdm a potential dis- 
tribution ({) as (ref. 20, 21) 


t". • . .. - -V?'. 

Th^/steady -velocity sat^i^^^^kontihuity equatfFn/(l45 ; hence, 

substitfc^hg equatidh'^t44TWTd5,^:*’" 
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\ys*: 


(45) 


L((J>) = 




(p?i).= 0 • 


Therefore, the correction velocity field satisfies a Poisson equation that 
becomes homogeneous in the limit as.pVj becomes divergence-free, i.e.,. equal 
to pUj. The boundary eendi:ti on for solution of equation (45) is expressed in. V. 
the inner product of equatiba’'{44) with’”‘the’' outward pointing normal to the 
solution domain closure 9fi._Since v-j is required to satisfy the boundary 
conditions for u-j , this yields 




(46) 


Therefore, equations (45)-(46) define the Newmann problem; <j) must be set to a 
constant at least at one nodal location on 9^2 to render the solution unique. 


Since the differential specifications for particular pressure and pertur- 
bation velocity are elliptic, their solutions are obtainable upon convergence 
of the parabolic equation solutions at any (each) axial station. A dominant 
term in the non-homogeniety of equation (45) occurs for i = 1, which corresponds 
to the axial derivation of the u^ velocity component. This is evaluated at the 
current Xi station using a second-order accurate backwards difference formula. 
Upon determination of the distribution of (j), the particular pressure equation 
is solvable since the divergence-free field ptj.j is known. At this point all 
dependent variable distributions are known at the current station. Based upon 
the previous station data, distributions for all variables at the next Xi 
station are predicted using a second-order accurate difference formula. The 
predicted values of parabolic variables are corrected by solving the initial - 
valued differential equations. Following convergence, the remaining variables 
are then corrected by solution of the elliptic differential equations. 


Coordinate Transformation For Juncture Region Geometry 

At this point in the development, the juncture region is assumed formed 
by two non-lifting intersecting surfaces with coincident leading edges, 
denoted Sj and $2 in Figure 5. The viscous juncture flow boundary-value 
solution domain is bounded by the envelope of Sj and $ 2 . denoted f£i(xi), 

the viscous intersection with the inviscid freestream, denoted f£2(^l)» 

two straight lines connecting their extremities, as noted in Figure 5^ The 
surfaces Si and $2 are not coordinate surfaces of x-j ; however, the fj^ are 
simply displacements parallel to the appropriate coordinate Xj^. For the 
transformation n.j = nj (x-j, f^i) that normalizes coordinates is 
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Figure 5. Juncture Corner Formed By Intersecting Parabolic Arcs 



(47) 


Here, the are assumed (at least) piecewise continuous curves describ- 

ing the closure and are normalizing coefficients. For example, shown in 
Figure 6 is a two-dimensional representation, where f2i is piecewise 
continuous while f22 is smooth. Using the chain rule, differentiation on Xi 




[j’22 ^ 


3X2 ' 

-3— = h - 9 
8X3 "31 8113 



(48) 
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Longitudinal Coordinate 


^1 


Figure 6. Coordinate System Transformation 


The functions , 1 £ i £ 3, are related to the metric of the coordinate 
transformation, equation (47), and are defined as 







(49) 


The superscript prime denotes the ordinary derivative with respect to x^. 

The n. coordinate system is fixed in the transform space and hence insensitive 
to . This non-orthogonal coordinate transformation is generally valid for 
3DPNS solutions provided the mean flow velocity vector u-j remains essentially 
parallel to n^- A curvilinear orthogonal transformation would be required to 
describe flows over highly curved surfaces using 3DPNS. 


Boundary Condition Specifications 

Figure 7 illustrates the viscous flow solution domain closure on 
the X£ plane. Pressure coupling with the three-dimensional inviscid 
exterior flow is accomplished by equating the complementary pressure Pc to 
pi on closure segment d-e-f. Furthermore, p^ = pi on segments c-d and f-a 
since these are assumed sufficiently remote such that boundary layer flow 
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exists thereupon. The normal derivative of vanishes on segment a-b-c. 

The particular pressure Pp is zero on segment c-d-e-f-a, and its normal 
derivative is constrained^according to equation (42) on a-b-c. 

No-slip boundary conditions are applied on the juncture region surface, 
f^l, hence segment a-b-c in Figure 7. The X 3 -derivative of both Ui and Vg 
are assumed to vanish on c-d, and U 2 can be solved directly from the primitive 
form of the continuity equation (24) as 

^ (50) 

3X2 3Xj 

Similarly on segment f-a, the X2-derivative of Ui and V2 vanishes and 


3-(|v^) _.^-a^p,u.i) 

■fix.3 V./. _ 


(51) 


Note the righft-band side of ‘botH* 'equations (50)-(51) contains the axial 
derivative of 'pu(, 'which' is formed .qs.ing a second-order accurate backwards 
•• difference formula. The continuity equation and the irrotational constraint 
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for the inviscid exterior flow are used to derive gradient constraint boundary 
conditions for Vjj, on the freestream closure segment. On e-f, this yields 


3V2 ^ 9V3 
9X3 9X2 


(52) 


3(pV 3) ^ _ 9(pUi) _ 3(pV2) 

3X3 9 Xj 9X2 


(53) 


On segment d-e the 2 and 3 subscripts in equations (52)-(53) are permitted. 
The normal derivative of the perturbation velocity potential function (ji must 
vanish everywhere on 9R^, as described previously except for being set to 
zero at one location. 


The boundary conditions for turbulence kinetic energy and dissipation 
function are vanishing normal derivative on segemtns c-d and f-a, and identical 
vanishing on the freestream segment d-e-f and the wall a-b-c. An adequate 
resolution of the transitional sublayer appears computationally uneconomical 
in terms of required discretization. Therefore, mixing length theory and 
definition of a dissipation length scale are employed to establish k and e 
at computational nodes lying closer to the wall than a specified level of the 
turbulence Reynolds number y+, 

pu^x 

y+ = -I - - - (54) 

U 


In equation (54), x^ is the coordinate x^ 
is the friction velocity 


T 


U. p 


perpendicular to the wall, and u^- 

(55) 


where is the Ui wall shear stress. In the near wall region, mixing length 
theory (ref. 15) yields the alternative expression for 


eff 

V 







(56) 


where to is Van Driest wall damping and a is the mixing length. Comparing 
equations (56) and (32) yields, for x^ £ y+ 
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(57) 


r ~ 202 '^^^ 

C 4 — ~ 


9X„ 


The additional required relation is definition of the dissipation length 
scale i, 

JIh = C V- 

d V e 


(58) 


which, for j< is assumed of the form 


KX 

V n 


(59) 


In equation (59), k is the von Karman constant (0.435). 


Finite Element Solution Algorithm 


The 3DPNS partial differential equation system for the viscous juncture 
region flow is identified. Each member of the equation set (25)-(29), (39), 
(43) and (45) belongs to the general class of second^order non-linear elliptic 
partial differential equations. Hence, for {q> h {uj, v^, k, e, p^, p , <))} 
each equation is of the form, ^ 



(60) 


In equation (60), the tensor indices range 2 t £ 3, and 1 £ i £ 3, K is the 
diffusion coefficient, fi is a function of its argument that specifically 
includes three-dimensional convection, p is any solution parameter including 
another dependent variable, and f 2 is the initial -value operator if present. 
The three-dimensional solution is required established on the bounded open 
domain fi 


fi = R^x 


Xq.X )e X 


Xi(0),Xi) 


(61) 


Note in equations (39), (43) and (45) f 2 vanishes identically and x^ appears 
only as a parameter. The general functional constraint for all q on the 
closure 9fi of fi, i.e., 9fi = 9R^x[x , x) e x^ x xi is 
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^-(p) * i'j q + 32 + a3'»*0 (^62) 

l 

where the a^. (Xj^,Xi) are specified. Except for q = {p(~, Pp, <|)}, an initial 
condition is required on x x^ (0) as 

q(Xi(0), X 2 , > 3 ) = qr,(X2«X3) (63) 


The finite element algorithm assumes all q and p interpolated on as 

q^Cx^.x) » (64) 

Determination of the expansion coefficients is accomplished using the Method 
of Weighted Residuals (ref. 22 ) as 


{N(x^)} L(qe)dT"- X 


{N(Xj)}ji(q_)dT 




|R| 


Rin9R2 

e 


im 


(65) 


where Sg is the assembly operator. The assembled finite element algorithm for 
the representative partial differential equation system is then 
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The rank of the global equation (66) is identical with the total number 
of node points on at which the dependent variable requires solution. 

For f 2 non-vanishing, equation (66) is a system of first-order ordinary 
differential equations. For f 2 vanishing identically, it is large order 
algebraic, and the matrix structure is symmetric ^ sparse and banded. 

An implicit solution algorithm is employed to solve equation (66) for 
all initial -value problem descriptions. The explicit form for equation (66) 
is 


i 






S {0} 


(67) 


The superscript prime denotes the ordinary derivative with respect to the 
Xi (m) coordinate. The matrix expressions for the first three terms, which 
account for axial convection and transverse plane convection and diffusion 
respectively, are 

* [C]p = tpui>l tN} {N} {N}'''dT 





{N> {N} l^'dT 


( 68 ) 




(N) 


3{N} 

3X, 




3x: 


&t 


Summation on 2 ^ £ 3 is implied, and derivatives on Xj are replaced by 

equations (48) in the sequel. All terms not explicitly involving q, i.e., 
{Q>, are contained in {fig- 

The single-step implicit integration algorithm employed is 


• m, + hl8{Q};,T + (1 - 9){Q> 


'4 


j+i 


'4 


(69) 


In equation (69), j is the axial station index, h is integration step-size, 
and e equals one-half yields the trapezoidal rule. Following the usual 
matrix manipulations, insertion of equation (67) into (69) yields a large 
order, non-linear algebraic equation system. The Newton matrix iterative 
algorithm for solution of this system is employed as 
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(70) 


0 {Q} 


j-i 






The dependent variable in equation (70) is the iteration vector, and 


(Q)S:i i {Q)j^, * mf,] - (71) 

where p is the iteration index, and £ > 0 is an integer that retards 
evaluation of the Jacobian as an economy measure. The right side of equation 
(70) is the homogeneous form of equation (69) evaluated with the pl!l iterate, 
i . e . , 




[C]. 


- (Q)jl * h 


* (1 - e)r|^>j| 


(72) 


where 




P = 


(QJ? + ff). 


(73) 


Equations (72)-(73) are defined only in terms of inner products on elements, 
with the assembly operator yielding the equivalent global expression. The 
vanishing of {F}, to within definition of a computed zero, yields equation 
(70) homogeneous, hence convergence of the iteration for any evaluation of 
the Jacobian. The initial estimate (Qli.i for any iteration can be deter- 
mined using 6 = 0. 

By definition, the Jacobian [J] for equation (70) is the derivative of 
equation (72) with respect to Each of the finite element matrices 

[C]e> [U]g and [Kj^ have imbedded therein functional dependence on each 
{Q}P^ 1 , see equations (68), and an exact evaluation of [j] is difficult. 
However, the exact Jacobian is not required by the formulation; an economical 
approximation useful for all five dependent variables is (ref. 23). 


[J] ^ 





- Me) 


(74) 
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The [C]g and [11]^ element matrices contain the exact non-linearity due to 
convection, while [K]^ accounts for the essential diffusive non-linearity as 
approximated by the effective viscosity. Dependent upon the overall behavior 
of the solution of equation (70), equation (74) can be updated every iteration 
with evaluated once every step with {Q}j> the previous converged 

iterate, or maintained at the evaluation of any converged iterate, 

In any case, the operations in equation (74) are limited to matrix inner 
products on an elemental basis. The rank of [Jj equals the order of {6Q}; • 
specific (Dirichlet) boundary constraints are applied within the evaluation 
of {F}. 
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NUMERICAL RESULTS 


Solution Initiation 

The basic requirement is to assess accuracy and stability of the devel- 
oped code for the 3DPNS solution algorithm. The test case corresponds to 
onset flow at zero angle of attack into a right juncture corner formed by 
coincident leading edge surfaces. The viscous displacement surface for tiie 
exterior three-dimensional potential flow is assumed the right intersection 
of two 10% thick parabolic arcs. The Hess computer program (ref. 14) was 
executed at NASA Langley using a (fine) discretization containing approxi- 
mately 5000 surface panels to produce the inviscid pressure distribution. 
Figure 8 shows computed spanwise distribution of potential Cp at three chord 
stations. The influence of the corner at the first two stations extends about 
0.3 y/C spanwise. The 3DPNS evaluation was conducted using the inviscid Cn 
distributions within the interval 0.04 < Xi/C < 0.14, wherein influence 
of the corner juncture was dominant within the first quarter chord. 

Figure 9 illustrates a representative 3DPNS corner discretization of R^, 
plotted without diagonals. The first step in the solution is to determine a 
sufficient number of complementary pressure distributions, since these depend 

only upon C (x , Xi), see equations (39)-(40). Core storage is minimized 
P ^ 

by solving only three distributions of p^(x„, x^) and interpolating to eval- 

C Jo 1 

uate both p_(x ) and the axial pressure gradient ap /ax^ distributions at Xi 
stations of the 3DPNS solution. Figures 10 and 11 illustrate the computed 
(dynamic) p^ and xi-derivative of p^ distributions at the initial station for 
the test case. Each row of data corresponds to the nodal distribution on 
"columns," extending from the wall to freestream, starting at the top-left 
extremity of Figure 9. 

The 3DPNS solution is initialized at a selected Xi(0) by using Cole's 
law (ref. 15) to imbed a turbulent u^ velocity profile on each node "column" 
with freestream value matched to the inviscid Cp. Figure 12 illustrates the 

initial Ui so obtained for the test case. Initial distributions for k and e 

are self-generated by the program using the initial Ui distribution, by 

extending onto R^ the mixing length concepts employed for the transition layer 

model. The sole requirement is to expand the dissipation length scale defi- 
nition, equation (59) to the local boundary layer thickness 6. This is 
accomplished in the program using an exponential decay term such that 
^c|('^) ~ equations (58)-(59) are combined with equation (57) to 

yield a solution for k and e. The initial distributions of k and e for the 
test case are shown in Figures 13-14. 

The transverse velocity field plus the <() and Pp solutions, cannot be 
evaluated at the initial solution plane, since each requires evaluation of 
9ui/8xi. Therefore, u and p must be assumed initially zero. 

X, P 
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Figure 11. Computed Complementary Pressure x. Gradient Distribution 
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Figure 12. Initialized Distribution For u. From Cole's Law 



Figure 14. Computed Initial Distribution Of Isotropic Dissipation 
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The 3DPNS solution is initialized with Uj.k, e, and four small (10~*^) 

integration steps are taken forward to generate sufficient Ui data to 
evaluate axial derivatives using a second-order accurate backwards difference 
formula. At this station, an estimation of v is computed by direct inte- 

gration of the continuity equation (24) in the form 

9(pVn) 9(pUi) 

” ' axi (75) 

In equation (75), x^ is the coordinate parallel to a node "column," and 
related to x by the familiar direction cosines a , i.e. 

aj MX 


^n ^n£^£ 2 < £ < 3 


(76) 


The solution to equation (75) is the scalar component of pVj^ parallel to 
the node column; hence, the Cartesian components are determined as 



V 

£n n 


(77) 


The V distributions are computed as this for the next four stations, 

X 

and are inserted non-iteratively in 3DPNS solution for Ui, k and e. Then, 
at stations 8-10, the transverse velocity solutions from equation (75) are 
iterated directly with Ui, such that a divergence-free three-dimensional 
velocity field pu^. becomes established. Hence, the transverse momentum 

equations for v^ are initialized with u^ at station 11, and the solution of 

the complete 3DPNS system initialized. Figures 15-16 show the computed v 
distributions at station 11. 


This completes the program generated solution initiation sequence which 
is descripted in program input terms in the following section. The last 
remaining requirement is specification of the turbulence model constants. 

The various coefficients are all NAMELIST input; their specification for the 
test case are, C<j)i = 2.8, C<j )2 = 0.45, C, = 1.0, C = 1.3, = 1.44, = 1.92, 

and C = 0.09. k e e e 
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Test Case Results 


Upon completion of the 3DPNS solution initiation, a specified number 

(four) of integration steps are taken at ax = 10 , to homogenize the v 

3DPNS solutions and initialize the predictor-corrector sequence for 
perturbation velocities and particular pressure. Figures 17-18 show the 
resultant computed perturbation velocity distributions, while Figure 19 
presents the particular pressure distribution at Xi/C = 0.046. 

Figure 20 is a victor plot of the transverse velocity distribution at 
this station. Note that the initiation procedure induces no tendency to roll 
up the transverse field. Figure 21 shows the corresponding computed particu- 
lar pressure, and dynamic complementary pressure distributions parallel to the 
corner bisector and parallel to one wall. The corresponding transverse plane 
gradients of particular pressure are almost an order of magnitude larger than 
those associated with the dynamic complementary pressure. However, absolute 
level of particular pressure is only ~ 2% of the dynamic complementary 
pressure. 

Figure 22 is the vector plot of transverse velocity distribution computed 
at Xi/C = 0.12, the station at which the 3DPNS computation was stopped. The 
combined action of the Reynolds stress distribution and the transverse plane 
pressure gradient has induced a reversed flow region emanating from the corner. 
A modest rotational region separates this reversed flow from the axial pressure 
gradient-dominated mass influx into the corner. The location of this region 
proceeds in a symmetric fashion away from the corner and parallel to the main 
diagonal. There is also an induced flow parallel to both walls and away from 
the corner. The net action on the predominant velocity Ui is to cause an 
undulation in the profile parallel to the wall, with the minima at the corner 
and at the lateral outflow boundaries. Figure 23 is the output from COMOC III 
for the computed ui distribution, and Figure 24 is the corresponding particular 
pressure distribution. 

Additional output from the code includes all other velocity and pressure 
fields, the turbulence kinetic energy and dissipation distributions, as well 
as integral evaluation of skin friction and displacement and momentum thick- 
nesses. These boundary layer integral parameters are evaluated by integrating 
parallel to node "columns" of the discretization. Upon a restart specifica- 
tion, the computation of the six components of Reynolds stress -u .u . using 

equations (30) can be commanded. Figure 25 is a plot of the distributions of 
uTuT at Xi/C = 0.12 along the mid-discretization node-row parallel to the 

wall. This is the transverse location of the extremum boundary layer shear 
stresses ufu^ and ufus on the farfield outflow boundaries. Both are character- 
ized by relatively flat plateaus with a sharp decrease of two orde rs of magni- 
tude about the corner bisector. The transverse plane shear stress u^ug is up 
to ten times smaller than these and exhibits a smooth increase to an extremum 
at the domain diagonal. The three normal stresses are nominally uniform with 
the Xx-component about twice the level of the transverse plane components. 
Figure 26 is a reproduction of the COMOC III computation of the complete dis- 
tribution of u^uT^ at this station. 
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Figure 21. Computed Initial Particular and Complementary Dynamic 
Pressure Distributions on Transverse Plane 
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Figure 26. Computed Reynolds Stress Distributions, 
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Figure 26. Computed Reynolds Stress Distributions (Contd) 
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Figure 26. Computed Reynolds Stress Distributions (Conld) 



DATA DECK SPECIFICATIONS 


Input facilities for the COMOC III computer program are sophisticated and 
greatly simplify data deck preparation and modification. The program sequen- 
tially scans the data deck and operates on command data cards as they are 
encountered. Numerical data required for each command operation are input in 
free format on cards directly following the command card. Command operations 
can cause arrays to be filled, initiate a series of solution operations or 
specify output formats and titles. Command card sequence is quite flexible 
and care has been taken to ensure that most operations which must be performed 
sequentially are specifiable under one command name. 

Most numerical data may be input in free format. Data delimiters may be 
blanks or commas, thus allowing for esthetic and meaningful arrangement of data 
and simple addition or deletion of interspersed numbers. Several features of 
free format input which greatly simplify repetitive and sequential data speci- 
fications are: 


Repetitive Numbers : 
Fills Array 

Repetitive Sequence 
(one per card only) 
Fills Array 


12.5*7 

12. 7. 7. 7. 7. 7. 
2(5. 2. 4. 

5. 2. 4. 5. 2. 4. 


Skip P locations 


10. 12. 3*P 22. T 


Fills Array 


10. 12. V V V 22. 



Increment by a constant 
Fills Array 

Exponential Notation 
Fills Array 


5*50 10 T 
10 60 110 160 210 

6. 10.0 E-2 14.E-4 T 

6. .1 .0014 


The data deck for a wing-fuselage juncture region solution is segmented 
into several sections for description, but appears in sequential order as 
indicated by the line numbers. The data has been lumped into meaningful 
categories, each of which can be related to the juncture flow problem and/or 
specific requirements for differential equation solution (eg., initial 
conditions, boundary conditions, etc.) The Fortran MAIN program contains 
important dimensioning parameters together with coding which links 3D inviscid 
interaction data with the viscous solution and is listed following the data 
deck description. 


I NAMELIST (Integer) INPUT 


2 

??opn‘; 


4jrn«(0 

5 WTNG-FUSPLAGE JUNCTURE FLOW. 


•_ ■ 


7 6N4MF01 


R 

NCOE » 210, 

LCOL a 50, 

KROW a 50, 


P 

NSNOOr a IP, 

NSEL=N a 2, 

TSIOE a A, 

NVAR a 3, 

10 

LG » 22, 

N5C a 0, 

NWALLS a 4, 

NTCHEK a 3, 

ll 

N«C = P, 

N'QKNN a 5, 

NEC400 a -A, 

NIZS a 250, 

1? 

NYY = 4, 

N2 7. a 4, 

JOROER a 1, 

NBAND a 23, 

i? 

NSpcr a u. 

NVP a 13, 

NSTRT a It, 

NOP a 10, 

lA 

NFt-^P = 1, 

TTK5 a 0, 

NOERIV a 2, 


15 

NO’-' = 101, 

NPTM a 200, 

NSTRFP a 1, 

NOOeS a 170, 

16 

NPVSX a 570, 

NOPVSX = 1250, 

NPVSXT a lA, 

IARR4Y(3P4l-3, 

17 

NVRHS a T, 


MLTLBS a 3, 

NVRH a 7, 

10 

TPTOOV a 1, 

tAPR«Y«28H=P048, 0, 50, 

62, 

IP 

NOP a 7, 

IPHI a t. 

NPRFS a -1, 

K8SAV =,-l. 

20 

I'‘CON a 1, 

L?HT a -1, 



21 

TRl a 1, 

N30PNS a 1 , 

NMCNTR a 1 , 

N2WAKE a 0, 

’2 

NMQU-^ a 2, 

NMROtr' a 50, 

NC a 8, 

NOUTS a 1, 

23 

Nri-MCG a 35, 

iriFRT = 30, 

NMOL a 8, 

NTMPLT a 1, 

26 

NPUTVC a JO, 

KNTP4S a 10, 



25 

26 S=NO 

KC'JMP a 1, 





Line Command 


3 


Solution Procedure (3D parabolic Navier-Stokes) 

5 


Data Deck Titles 

6 

FENAME 

Command to initiate Namelist read 

8 

NODE 

Slightly larger than the number of nodes in solution 

9 

NSNODE 

Number of nodes in MACRO element geometry 


NSELEM 

Number of macro-elements to be refined 

10 

LG 

Set greater than the total number of one-dimensional 
node sets specified under command CNTNDS. 


NBC 

Number of gradient boundary condition sets 

11 

NEQ 

Number of variable arrays in solution (LISTED following 
command IPINT) 


NEQKNN 

Number of variables being integrated (eg. first 5 
listed following IPINT) 


NEQADD 

Number of variables not being integrated during first 
few initializing integration steps. 

12 

NBAND 

Jacobian band width (must be an odd number) 
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I NAMELIST (Integer) INPUT (contd.) 

Line Command 

14 NE1E2 Turbulence Model Switch 

0 - No MLT 

1 - MLT 

2 - MLT when current station is greater than 

E1E2SW 

ITKE Turbulence Model Switch 


0 - No TKE (see C4EDSW in NAME02) 

1 - 2 eqn. TKE, DISS Model 


16 

NPVSX 

Number of entries in Cp tables (formed) 


NPVSXT 

Number of entries in Cp tables (inviscid data) 

17 

NVRHS 

Perturbation pressure variable No. (Parameter arrays 
must be in order; Pp, p^, 4>) 

20 

ICORN 

Set to 1 for juncture region flow 

21 

NMOUT 

Standard print format 



2 - Arranged in geometrical shape (some data may be 
omitted) 



3 - Columnar by node number 


NMBOUT 

Number of variables to be printed (specified under 
lOSAV) 


NC 

Number of significant figures printed 

23 

IDIFRT 

Debug from WLFLXS (number of times) 

24 

KNTPAS 

Causes print to occur as a function of the pass counter 
(over-rides DELP) 

25 

KDUMP 

ECHO print of input data 
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II NAMELIST (Real) INPUT 



r.N»«E02 








!JIM= = 300., 

TOFINP • 530., 

PINF • 2116.8, 

REFL - 

1.0, 

29 


rCN » .41. 

C"MPX * 2., 

COHPY « 2., 



30 


0»'!4*Ya26»»l. 

,1.0^6, 

3*1.0, 





CF = 1.0, 

CW - 1.0, 

YPLUS • 30., 



22 


YT = 2.0. 

5LTH = 0.5, 




33 


SIMOLT « .01, 

CHIF05 = l.e-5 

.CHISTP « 10., 

TMULT • 

' 1.44, 

34 


OS’^A’T . 1., 

«:r=25w= 10000., 

C4EDSW=.C0O99, 

VSTART 

» 0. , 

35 


Fi:miiL'-=-.002 , 

HSTNIT»l.0'-4, 




^6 


TO = .11, 

TO » o.ia. 

DELP « 0., 

HMftX - 

2., 

37 








Line Command 


28 

UINF 


TINF 


PINF 

29 

COMPX 


COMPY 

33 

SIMPLT 


CHI EPS 


CHISTP 

34 

DSTART 


E1E2SW 


C4EDSW 

35 

HSINIT 

36 

TO 


TD 


DELP 


HMAX 


Velocity non-dimensional izing factor 

Temperature non-dimensional izing factor 

Pressure non-dimensional izing factor 

Direction 3 compression factor for geometric form 
print (NMOUT = 2) 

Direction 2 compression factor for geometric form 
print CNMOUT = 2) 

Station at which to switch from explicit to implicit 
integration. 

Matrix iteration convergence interval. 

Implicit integration step limit 

Integration station where NE1E2 flag is changed (MLT) 

Integration station where ITKE flag is changed (TKE) 

Initial integration step size 

Initial integration station 

Total integration distance (TF = TO + TD) 

Integration print interval (% of TD) 

Maximum allowable step size 
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Ill DYNAMIC ARRAY DIMENSIONING AND DEPENDENT VARIABLE SPECIFICATION 


BBIPINT -1 

41 

47 


15623789 10 0, 
177568 4*0t 
10411 1 T 


Line Command 

38 FEDIMN 

39 IPINT 

40 


41 

42 


Dimension arrays to fit problem size 
Cards following specify dependent variables 
Variable 1 parimary flow velocity 

5 TKE 

6 Dissipation 

2 Transverse velocity (V 2 ) 

3 Transverse velocity (V^) 

7 Perturbation pressure (Pp) 

8 Complementary pressure (p^) 

9 Continuity equation potential ^ (cjii) 

List of variables which are extrapolated at each 
new implicit integration station. 

Variable number assignment 


Note: NEQKNN in NAMEOl specifies the number of integrated differential 

equations (i.e. the first NEQKNN Variables in Line 40 are integrated). 
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IV FLOWFIELD GEOMETRY 


6?L INK A 

9 




aatappay 

A50 

250 

T 


A50ST3D 

-1 




A5 



1 -2 T 


47LTNR3 

lA 




ASNf-TA 





49 



9 9’’ 


509=05 





51 



9 9 T 


52S^YOo 





53 



2*A T 


'ASFLGN 





55 



9 10 3 

1 12 13 

56 



1 3 lA 

15 19 16 

57f>roy4i? 

280 

290 I2A8 


5 0 


0. . 

.8 .01918 

0. .025 

59 


.2 

1 . 25.8 . 

01918 0. 

60 


0. 

0. .01918 

.05 0. 

61 


0. 

.0009 .01918 .1731 

62 


0. 

0. .758 0. 

0. -5 .5 0. 

63 



.58 .720 

0. .71 T 

6AorNP 





65P=aO 

5 

63 

-26 


66 



11 1 2, 

CM 

670=40 

5 

63 

-26 918 


68 



182 181 171, 171 172 


19 Z 
17 18 T 

.05 .025 .0 .1731 .1731 

1.25 1.25 1.25 1.25 

.025 .05 .025 0. .01918 

.1731 .08 .0009 .09 .0009 

. 783 0. . 720 . 58 . 783 0. 


T TURN L0M9R RIGHT DIAGONALS. j 

1 

182 T TORN UPPER LEFT OIAGQN/^S. 


Line Command 

43 LINK4 9 

45 PSIBD 

47 LINK2 14 

48 NETA 

50 NEPS 

52 STYPE 

54 SELCN 

57 DEPVAR 


Dimension arrays required for grid generation 
Input intended diagonal skew of generated elements 

1. ) Positive Macro-Element No.- right running 

2. ) Negative Macro-Element No.- left running 

Reads discretization data and generates triangular 
finite elements in each specified Macro-Element 

Number of divisions in local coordinate direction 
(n) 

Number of divisions in local coordinate direction 

(e) 

Type Macro-Element 

3 = Triangular 

4 = Quadrilateral 

Macro-Element nodal connection table 
Coordinate and initial values to be distributed 
(Variables X 3 , X 2 , and Ui) 
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IV FLOWFIELD GEOMETRY (Contd.) 


Line Command 

64 DONE 


End geometric data. Grid is refined to specifica- 
tion and generated grid and finite element connec- 
tion tables are stored for reference. X 3 side node 

data is a geometric progression ratio (x) where: 


15 17 14 




- Primary Flow Direction 
Xo,X- - Secondary Flow Plane 

u 6 
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V BOUNDARY NODAL ORDERING ARRAYS 


69CNTPTS 

70 

77 

73' "OPO 

74S.00 

75 

i 77 LTNK 3 

I 73 T»®o»Y 

7 °LT')K 1 


-1 

19*10 T 

-1 

190*11 1 


10*11 It 18*n0 20t 

4 T 

250 0 T 

4 T 


,v. .. 

9*1-10 189, . 17*1-10 171 ,T • ’'"•A 

• .. 

N0N-D1M5NSI0N. CONSTANTS IsiSEN) '' ' . , ^ 

IPHt IN GCnMPL ' ■ 

FINITE ELEMENT MATRICIES tGEOHFU •- ±_ 2 __ ^ 


Line 

Command 


69 

CNTPTS 

Number of nodes in each set 

71 

CNTNDS 

Node numbers in each set. Nodes are ordered in 
sets normal to boundaries of interest 

73 

IBORD 

Sequential numbering of boundary nodes in order 
parallel to boundaries of interest 

77 

DIMEN 

Form non-dimensional izing factors 

79 

GEOMFL 

Form Finite Element matrices. 


Note: 1) LG in NAMEOl must be equal to or larger than the total number of 

sets of data. 

2) KWFLXS in NAMEOl indicates number of sets of this data over which 
to compute skin friction. 

3) LMLT in NAMEOl is the number of sets of data over which MLT model 
is used to compute turbulent viscosity. 
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VI TITLES AND HEADINGS 


Line 

83 

86 

90 

105 

109 

115 

119 

122 


soo*=sroti*^ 

81 W^NG-FUSELAGE JUNCT'JRS flow. 


83C-WTITLE 

6A WTNG-FUSFtAGF JUNCTURE FIGH. 

.8 5CCNF 

95''»"SrRTPT 20« T OFSCRIPTIVS TITLE AT BFGINING OF OUTPUT. 

87 WtNG-FUS'=LAGE JUNCTURE FLOW. 

PP 

PFOONF 

PO"':SrPIDT 33? T lOPAR PARAMETER TITLES FOR OUTPUT, 

qj OFcqjjc.;-F FNGLISH-FT FMGLISH-IN m-K-S C-G-S 





T 

M 






ki A 

•M/S •••••••« 




qC^iUTTV 

.LBM/FT3.... 


.KG/M3 

• G/CC • • • • 


95 

TrMPP9A'"U®F.*.. 

.RANKINE.... 


.KELVIN 

* N* &••••• 


96 

'ENTHALPY 

.R’'U/L8M.... 

*M«Aa« •••••• 


•N*A***«« 


97 

=RO?.SP'C, F’^&T 

,RTU/L'’«-R. . 

•NaAaa****** 

.KJ/KG-X.... 

• N* A • • • • • 



VT 

.LRM/F+-S.. . 

• M*A»* •••••• 

.NT-S/M2. ... 

.POISE... 


99 

lOrAL PPESStiR = 

,PSR 


• NT/M2 

.TORR.. .. 


100 

L''FAL SOLirirN 

.MACH. NO. 

•*DPDXX*«* 

..ENERGY.. 

.CON. VAR. 

101 

NWGFPM H'S 

• • • H? 1 • « • • 


• • • G23« • • • 


• 

102 

H*S 


• • • 2* • • • 

• • • G33 • • • • 

• • • *GX • • • 

• 

103 

XI /L REF 

•DXl/LRFF. 

.EPSILON.. 

.OXIH/LREF REEL 

REYNOLDS 

NO 




105WPA®A 

105 

107 

1C8 

109*0NUMB ■ 

no 

lU 

112 

■113 , 

-114 

1 1 EOF SORT P” 
115 UV/UREF 
117PP ; 

118PON3®* . . 
IIRT'-SavS 


120 

1?1 

I’ETC^ULT 


1 


1 


203 


I 


1 


5*2, 2*2 162 164 153, 3*2 164 153, 3*2 170 174, 

3*2 165 2, 2 -175 3*2, 3*2 176 2, 3*2 177 178, 

2 2 159 158 157, 3*2 108 2, 5*2, 5*2, 5*2 T 

999, 5*200, 999, 200 4*43, 200 27 200 2+27, 

200 10 200 2*10, 2C0 58 200 58 200, 

200 97 200 97 200, 200 30 200 30 200, 

200 38 200 2*38, 999, 39 4**6, 300 154 100 135 122, 
200 A-*!! 186, 200 4*11 139, 11 12 14 85 47 T 
T IFVTHD ttTL'S for OUTPUT CEPENDENT VARIABLES, 


V2/1JPFF 

V3/URFF 

V2 PR 

V3PR 

PHIl 

TKE/TKERFF 

CISS/DISSREF* 

NUANUREF 


1248 

2248 3248 

260 

262 

7248 

9248 5248 

6248 

1247 T 


123 


9*2 21, 10*1 T 


124PLOTS -1 

125 2248 3248 1248 t U2, U3, U1 


Command 

COMTITLE 


MPARA 

lONUMB 

lOSAVE 

lOMULT 


Problem identifying titles 

Std print titles 

Std output header labels 

Scalar parameter print multipliers 

Scalar parameter print locations 

Dependent variable print labels 

Dependent variable print locations 

Dependent variable print multipliers (RAgRAY LOC,) 
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VII DEPENDENT VARIABLE AND PARAMETER BOUNDARY RELATIONS 


Line 

129 

133 

153 


Note: 


I?6T 

\ 27 ~ — 
17<)KR»)C 


PCUNPiRY rONOI'^TONS 
1 


’31 




19*110 

1 

T 

FIX WALL NODES. 

132CCNP 








133KRNT 

2 


1 

T GRADIENT B 

. C. CARO 0 0 0 AI 2 A3 2 • 

134&DC 

0 

0 

0 



0. 

2 0. 2 

135 




19*110 

10 

T 

B. C. ON FRFESTRFAH NODES. • 

136400 








137 




19*110 

1 

T 

FIX WALL NODES. 

13 800NS 








1 ’gKONO 

3 


1 

T GRAPIUNT B. 

. C. CARO 000 AI 2 A3 2 

140400 

0 

0 

0 



0 . 

2 0.2 

14t 




19*110 

10 

T 

8. C. CN FR'ESTREAM NODES. 

142400 








143 




19*110 

1 

T 

FIX WALL NODES. 

144f'0»)= 








145K0MC 

5 







1454'V’ 








147 




19*110 

1 

T 

FIX WALL NODES. 

i4acrN= 








140i<BM3 

6 







150400 








151 




19*110 

1 

1 

FIX WALL NODES. 

15 200M5 








153KR4I0 

7 


1 

T GRADIENT B, 

, C. CARO 000 AI 2 43 2 

154400 

0 

0 

0 



0« 

2 0.2 

155 




19*110 

1 

T 

B. C. ON WALL. 

1564C5 








157 




9*11 2, 

17*110 

20, 9*1-1 190 T FIX OUTER BOUNDARY. 

! 55'^'IV = 








1 50564''’ 

8 







1604''0 








IM 




9«I1 3, 

17*110 

20, 9*1-1 190 T FIX OUTER BOUNDARY. , 

le200N= 







, 

16‘’5'<'!0 

9 







1644C’ 




DON' 




165 




19*110 

1 

T 

FIX WALL NODES. ^ 

Command 








KBNO 1 




Nodes 

where 

Ui is held constant at initial values 





(Wall 

Nodes) 



KBNO 2 


KBNO 7 


U 2 normal gradient boundary conditions (aa inter- 
nally generated). Followed by nodes where U 2 is 
held constant. 

Perturbation pressure boundary values (internally 
generated) 

1) ai = 0.0 indicates symmetry plane boundary 

2) ai = 1.0 indicates no-slip boundary 

3) ai = 2.0 indicates slip plane boundary 


1) NBC in NAMEOl is set to maximum number of boundary nodes. 
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VIII TABLE DATA 




laTCN" 

16<?T 

17C- 

X 

Y 

30 INVISCIO SOLUTION ( 
2 

FIRST ITERATION ) 

CP 


171 

. 1C7600 

0.0 

.010180 

.999000 

.077580 

.156600 

-.03 



.107600 

.027840 

.019180 

.999000 

.077580 

. 156600 

-.028630 

625 

171 

.1C760Q 

. C47310 

.019190 

.'991300 

•113600 

.155400 

-.019750 

626 

174 

.107600 

.0 71650 

.019180 

.991800 

.107100 

.155400 

-.019220 

627 

175 

.107600 

. 102100 

.019180 

.995000 

.093740 

.155900 

-.023190 

528 

1 76 

.107600 

. 140100 

.019183 

.999600 

.073910 

.155700 

-.030030 

629 

177 

.107600 

.190000 

.019180 

1.004000 

.063880 

.157400 

- .03 785 0 

630 

178 

.107600 

. 268700 

.0191 «0 

1.009000 

.047310 

.158200 

-.045550 

531 

179 

.107600 

.418400 

.019180 

1.C12000 

. 028460 

. 158600 

-.050210 

632 

180 

,107600 

. 767700 

.019180 

1.013000 

.010420 

.158200 

-.044320 

633 

181 

.107600 

1.5160.00 

.019180 

1.004000 

.0021 83 

.157400 

-.033560 

634 

181 

.107600 

2.013000 

.019183 

1.301000 

.000353 

.157000 

-.027620 

635 

18? 

.107600 

6. 008.000 

.019180 

. 998600 

.000074 

.156500 

-.021580 

636 

184 

.107600 

9.002000 

.019183 

.998500 

.001582 

. 156500 

-.021550 

637 

185 

.190700 

0.0 

.030820 

1. 108C06 

.061740 

.137200 

-.265000 


1 86 ; 

.190700 

.039460 

.03082P 

l.icaoco 

.051740 

.137200 

-.249200 

654 

187 

.190700 

. 058910 

.030870 

1.C99000 

.100900 

.136100 

-.235900 

56 5 

188 

.190700 

.08’220 

,O-’08 7Q 

1.094000 

.098580 

.135500 

-.225100 

666 

139 

.193700 

.113600 

.0208 70 

1.093000 

.089010 

.135100 

-.214700 

557 

1°0 

.1907C0 

. 151600 

.030870 

1.087000 

.076350 

.134700 

-.205400 

668 

191 

.196i700 

.201400 

.030820 

1.CE4000 

. 062100 

.134300 

-.197100 

669 

192 

. 19C7C0 

.280000 

.030870 

1.081000 

.045330 

.133900 

-.187900 

670 

19? 

.190790 

.429600 

.030820 

1.075000 

.026180 

.133300 

-.175100 

671 

194 

.190700 

.778500 

.030870 

1.067000 

.008928 

.132200 

-.156500 

672 

195 

.190700 

1. 52 6000 

.030870 

1.050000 

.001769 

.131300 

-.140500 

573 

196 

.1 90700 

3.022000 

.030820 

1.055000 

.000281 

. 130900 

-.133100 

6 74 

197 

. 19C7C0 

6. C12000 

.030820 

1.054000 

.000059 

.130600 

-.128200 

675 

198 

.190790 

9.C03000 

.030820 

1.C53000 

.001325 

.130500 

-.125900 

676 

199 

.290800 

. 0 

.041193 

1.191000 

.041440 

.099540 

-.442 

703 

700 

.290800 

. 049830 

.0411°0 

1.191000 

.041440 

.099540 

-4429500 

703 

201 

.290800 

.069250 

.041 190 

1.1E3C00 

.073830 

. 098910 

-.415500 

704 

202 

.790800 

.093540 

.0411'50 

1.176000 

.073420 

.098310 

-.398200 

705 

70? 

.290800 

-123900 

.041193 

1.168000 

.C6745G 

.097630 

-.378300 

705 

204 

.290800 

.161800 

.0'Vn93 

1.159000 

.058860 

.096910 

-.357100 

707 

705 

.790803 

. 211600 

.041193 

1.1 50000 

.048500 

.096140 

-.334500 

708 

?06 

.290800 

.293200 

.341190 

1.139000 

.035500 

.095220 

-.3C8000 

709 

2C7 

.290800 

.439500 

.041190 

1.126000 

.020150 

.094080 

-.276000 

710 

208 

.290330 

.788100 

.041190 

1. 11 0000 

.0C65C4 

.092810 

-.241500 

711 

209 

.290800 

1.535000 

.041190 

1.101000 

.001227 

.092000 

-.219800 

712 

710 

.290800 

3.029000 

.041190 

1.09 70-30 

.000192 

.091690 

-.211600 

713 

21 1 

.29CECC 

6. C160C0 

.041190 

1.095000 

.000040 

.091510 

-.206800 

714 

212 .290800 0.004000 

713VU2P05 -3 

’1 4 .' 

215 .i 

216 

217VU2V4L -3 

718 
2’9 
220 

22n"J3POS -3 

22? 

223 

224 

725VI3V4U -3 

776 

227 

278 

.041100 

00999 .0194 
08459 .1076 
2559 .2908 

173 .17473 
1°41 .198 . 
2173 .2191 

00999 .0194 
CB459 .1076 
2559 .2903 

173 .17478 
1941 .198 . 
2173 .2191 

1.093000 .000948 .091380 -.203400 

ViP. GPr-M TA3LP STATIONS JXll 
.03126 .04628 ,06409 
.1331 .1609 .1907 .2224 
.3271 .3644 .4025 .4413 .4804 T 
^ VAR. GPOM. X2 TABLE 
.177075 .1798 .183 .1865 .1902 
2018 .2055 .2090 .2122 .2150 
.2203 .22088 T 

T VAR. G53M, TABLE STATIONS (XI) 
.03126 .04628 .06409 
.1331 .1609 .ljMA7 .2224 
.2271 .3644 .4C25 .4413 .4804 T 
VAR. GEON. X3 table 
.177075 .1798 .183 .1865 .1902 
2018 .2055 .2090 .2122 .2150 
.2203 .22088 T 

715 


Line Command 


167 END 

121-212 

213 VU2P0S 


Returns to Main and Cp table data is read. 

3D Inviscid Cp solution and coordinates 

Coordinates in the primary flow direction where X 2 
variation is specified (non-dimensional i zed by 
ALC) 
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VIII 

TABLE DATA (Contd.) 



Line 

Command 



217 

VU2VAL 

fz at each Xi coordinate followed by fi at each xi 
coordinate 

221 

VU3P0S 

Same as for VU2P0S (X3 

direction) 

225 

VU3VAL 

Same as for VU2VAL (X3 

direction) 


IX SOLUTION INITIATION 



Line 

Command 


229 

LINK2 30 

Sets up Pc tables 

230 

LINKl 11 

Evaluates p^ from tables at initial station TO 

232 

LINK2 10 

Generates Uj initial distribution from Cole's 
and Pc array. 

233-240 

LINK 

Initialize other parameters 

241 

LINKCALL 

Parameter evaluation at each integration step. 


1 11 Evaluate this step 

2 4 Evaluate second-order du/dx 

2 15 Evaluate flow integral parameters 

2 3 Evaluate wall shear stresses 

5 6 Evaluate turbulent diffusion coefficients 

5 5 Evaluate perturbation pressure (pp) 

243 QKNINT Initiate integration of dependent variables 

244 EXIT STOP 
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X FORTRAN MAIN PROGRAM 


0001 

0002 

W>(?^ 
SO 04 

0005 

0006 

0007 

300P 

0005 

0010 

noil 


C - 0 - M 


C 


COMMON tTSTtZOOl 

P0UIVSL5NCE ( LISTtOOtsn. JPR 


1 


COMMON / VJPBLr / IA1P, 4Y! 005001 . RAPRAY(00500) 

0IM5NS10N RZIlIt 1.(400) 

C0UIVAL5NC5 ( 5 7(11* 17(11* | iifn7m |7 1 

J^CUIVALPNC' (1(0711, IX > • <l-( 072 1 ,I Y 1 ’ lM?In ’ ^^pvsxl 

* ,(L(074).TC5 1 , (1. ( 139 I f IX3ST 1 1 (L ( 140 1 * ICPVSX) 

^OUI VALCNC*^ ( I A P 7 A V ( OOO 61 1 * KO'JMp 1 
roUtVALFN':'^ ( T4»P AY (000'>21, 17SI7F 1 
/ A°5AYS / 77(690001 
MZ5175 = 69000 

CALL =7PS"T ( 207, 500. -1- 1 - o. 5i 


•'. 77 '* 


0029 

0030 

0031 
0332 
C033 

0034 

0035 
0035 
0037 
00’ 9 

0039 

0040 

0041 
3042 

0343 

0344 

0045 

0046 

0047 

0048 


Line (ISN) 
3 
9 

10 

13-16 

19 

20-44 

45 

46 


200 

50C0 


300 


310 


(ICP + IMll 


-- nr*'*’* T I-- 1 1 

07 ((X 3 ST+ 7 _i ) = Pznx + I-17NPTM11 
cpTMlTI 3'7 10«6 *30X *910. 6 1 

T - 0 

30 300 K=1.NPT 
no ?oo j=i,3 
7 = T*1 

Lnp = NP^^U-l) + K - 1 
9 7 (T 0 0(79 X* 7- 1 1 = PZ(1C9+L0C1 
oo 310 J=1 fX7 
LOC * J 7 XT 

P7 ( t CPVSX + LOC-l 1 = RZHY+J-ll 
C0MT7M-J5 
NW = 6 

HO 7TO(NW, 61101 (07(lX3ST7l-l1*I’l*NX1 
H0TT9(NH,6n01 (R7( TCPVSX+XT + I-11 *1=1* XT 1 

vi07TP(MM. 61101 (97(7CPVSX + I-11.I*1*X.71 _ 

6110 POPMAT ( IH , 45H XI STATIONS, X2 VALUES, PRESSURE TABLE. 
1 (3E15.51 1 

CALL 531 NOT 
30 TO 100 
FNO 


Comman block for scalar arrays (Integer, Real) 

Common block for variable arrays 

Set program arrays size indicated equal to IZ array dimension 
Initialize all arrays to zero 

Begin reading data deck (return upon encountering command END) 

Read and store inviscid table data 

Continue reading data deck commands 

Begin new problem (terminates with command EXIT) 



SUMMARY 


A parabolized three-dimensional partial differential equation system has 
been derived for prediction of turbulent subsonic aerodynamic flows in junc- 
ture regions. A Reynolds stress constitutive equation is applied which intro- 
duces source terms into the transverse plane momentum equations. These terms 
depend primarily upon the shear of the axial velocity, and the levels of tur- 
bulent kinetic energy and isotropic dissipation. They exert a dominant influ- 
ence on the momentum equation solutions. The developed 3DPNS equation system 
and implicit solution algorithm have been coded and evaluated using the COMOC 
III computer program, and satisfactory performance has been achieved. The 
3DPNS equation system, upon pressure coupling within an iterative interaction 
solution with a three-dimensional potential flow code, can provide a computa- 
tional methodology for determination of juncture region geometry influences on 
flow structure and induced drag. 
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APPENDIX 


The 3DPNS option in COMOC III is implemented using linear interpolation 
polynomials spanning triangular-shaped planar finite element domains 
possessing vertex node points. The basic geometry is illustrated in Figure 
A.l including the various coordinate systems employed to span and R|. The 
normalized, 1 inearly-dependent coordinate system ?-,• is the fundamental^ 
descriptor for R| spanned by linear functions; hence (ref. 22) 


and k = 1 in equation (64). 



Figure A.l Coordinate Systems for Two-Dimensional Finite Element 
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The c and n coordinates are related as 


{c> M 


[€^in/|j(2)]c2/?2^3) 


Ci/q(2) 


where ^-j(n) is the i^h coordinate pair of node n. Furthermore, 


(A. 2) 




(A. 3) 


where the elements of a^-j are the direction cosines of the coordinate 

transformation that aligns with the line connecting the first two nodes 

of R2. 
e 

The gradient vector operation is fundamental to the formulation of the 
equivalent matrix structure. Using the chain rule and the summation 
convention. 


3^4 »n 




(A.4) 


The coordinate transformations, equations (A.4) and (48) facilitate the final 
two evaluations. The gradient operator on 5 -j is an elementary analytical 
operation using equation (A. 2). Derivatives on x^, see equation (48), **. 

introduce additional contributions. ; 

i 


The second calculus operation intrinsic to the algorithm is integration'* 
'Of products of the elements of {Ni(x )} over R^. For example, in the primary 

mean momentum equation (25), the i =1 convection term corresponds to f£ in 
the general description, equation (^^/ identified with u^. Assuming 
,fof simplicity that variable geomet»^,‘^sfilaiis’^only the X 2 coordinate, i.e., 

VfgV “ ^0^ ■'n equation (49) this term Keccfeies 


pu 


3Uj 


.^4 








®n2 





(A. 5) 
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Inserting equation (A. 5) into the algorithm formulation, equation (65), and 
„ using the elemental interpolation for both pui and Uj on and upon 

. i rearrangement of scalars, the algorithm equivalent of equation (A. 5) is,''J>*-- 


ii 


■ 


* 


■>zil 





Un) {rhoud 


; - ( 


^22 




^{N}>H0W^ 








(A. 


•' Recall that only the elements of {Ny}'are functions of n» and that the 
derivative of these linear functions* are constants on R^, 

Generalizing to 1 <. j :< 3, fornfation of the finite element equivalent 
of the three-dimensional convection term in equations (25)- (29) requires 
directional derivatives. Hence, equation (A. 4) becomes 


* i*2j(xj){Sll2}ge2 * ^3|(*|Htn3Je«3 


(A. 7) 


where the (fixed) coordinate transformation, equation (A. 3) is imbedded into 
the {BllJ}. J is the index synonymous with e-, the unit vectors of x^, and 
the elements of {BllJ} are independent of Xi.'^ Using these concepts and 
notation, the linear interpolation finite element equivalent of the entir^.o 
three-dimensional convection term imbedded' in each of the 3DPNS equation^/, 
-as represented within fn of the general ’differential equation (60) becomes' 




R2 

e 


{N)p5, H;dx = A ■! 


{RHOUDJ 83000] 




[8200] {RH0U2}g - g22{RH0Ul}g h2j{B112}g{Q}^ 
[8200] {RH0U3}g.,- g32{RH0Ul}g h3^{B113}g{Q}^ 


{RHOUDg [83000] g23{ETA2}gh2^{B112}g{Q}g 

J V, 

+ g33{ETA3}gh3^{B113}^ {Q}^ 


(A. 8) 


The elements of {ETAJ} are corresponding nodal coordinates of R^. Matrices 
with second index >2 are hypermatrices, elements of which are 
themselves matrices, see reference 22, and inner products must be performed 
in the correct sequence. 

The finite element form for the equation determining particular pressure 
distribution warrants comment. From equation (43), and enforcing the con- 
tinuity equation (24) directly, imbedding the resultant statement within the 
finite element algorithm, equation (65), and performing select integrations 
by parts, yields the following statement for particular pressure p (x ;x^). 


fa{N} ap. 




dx + 


apu. 8u. 

{N> dx 


R2 

e 






( 3(N} 5 


3Xj 




3(pu,) 

^ a + — L 

3X„ £ 3X, 


dx 


{N} 


au, 


pu 


j ax. 


ax, 


.au, 


au, 


ax. 


+ 


ax. 


3R2o3R2 
e e 


n^dx 


' 

o{N} 


1 

u -n • dx 


3X^ £ 3Xj^ 

:J J 

J 

O 

- 



E {0} 


(A. 9) 
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l^ation (A. 9) employs the limited summation convention, 1 £ i, j £3, and 
2 1 k, £ £ 3, and equations (A. 4) and (47) must be employed in the matrix 
development to replace derivatives on x^. Note in particular that 3{Nj}/axj 
does not vajiish identically in the variable .iuncture geometry. 






The first term in equation CA.9) is the Laplacian operator written on 
particular pressure. The second and third terms result from an integration 
Ipy parts of the second term in the parent equation (43), and simplification 
' ng the continuity equation (24). The fourth term is the boundary condit 
dement, equation (42), and valid on symmetry plane, no-slip and slip-wal 
re segments. The l-ast t^erm is a closed surface integral over an xi 
ation step, resulting frc^^use af a Green-Gauss theorem in the 
gration by parts. Since vanishes on all closure segments 9R| for 

confined flow, this term involves an integration over R| at the current 
station, followed by subtraction of the previous station evaluation on an 
assembled basis. 

Since the convection and boundary condition terms in equation CA.9} 
involve double summation index sets, the resultant expansions in transformed 

4 pordinates are very lengthy. Nevertheless, upon application of the finite _ 
l^ent interpolation, equation (64), employing the developed standard matrix 



forms., and combining like terms, each of the integral expressions can he 
djrectly evaluated on a triangular non-regular discretization of R^. For 
■ Example, for coordinate stretching spanning X 2 only, and recalling that 
< 0 < 3 includes initial -value derivatives on Xx, the second term in 
I’^'^^quation (A. 9) (9.9) becomes, 




3pu. 3U-* 

- * 


Cl 



{RH0Ul}e{B112}gh2i[[B200]{U2}g - g22<B10} - g23[B200] {ETA2} J^^CB wJ{U2}^ 
{R^l>e|B113}g[[b200]{U3}g - g22{B10} - g23[B200] {ETA2}gjh|j{ilrfig{U3}g 


+ {Bl^{RH0U2}^{B112}^h|^{iimJ{U2>^ + pOU3)^{B113}J{BU3>^{U3| 


m- 


+ 2[{RH0U2>e + {U2} J {B113|h2i^B112}j[{U3}^ + {RH0U3}^]J 




(A. 10) 


C^.IO) is a columri rfflitrix'of degree three, a,oxl. the matrices with‘«’» 
jrscr’ift prilfff Tle'note the nodal distfitiution of xi-derivatives of the 
Corresponding dependent variable." the coding of the PPRES subroutine in 
ulWOC III, for evaluation and solution of equation (A. 9), contains terms 
written in the formalism exhibited in equation (A. 10). 
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turbulent flow in the juncture of two intersecting parabolic arc 
airfoils. Specifications for data input requirements for COMOC III 
are discussed. 
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